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Abstract 



In this paper we present the results of a large-scale numerical in- 
vestigation of structural properties of a model of cell membrane, sim- 
ulated as a bilayer of flexible molecules in vacuum. The study was 
performed by carrying out extensive Molecular Dynamics simulations, 
in the (NVE) micro- canonical ensemble, of two systems of different 
sizes (2 X 32 and 2 x 256 molecules), over a fairly large set of tem- 
peratures and densities, using parallel platforms and more standard 
serial computers. Depending on the dimension of the system, the dy- 
namics was followed for physical times that go from few hundred of 
picoseconds for the largest system to 5-10 nanoseconds for the small- 
est one. We find that the bilayer remains stable even in the absence 
of water and neglecting Coulomb interactions in the whole range of 
temperatures and densities we have investigated. The extension of 
the region of physical parameters that we have explored has allowed 
us to study significant points in the phase diagram of the bilayer and 
to expose marked structural changes as density and temperature are 
varied, which are interpreted as the system passing from a crystal to 
a gel phase. 

1 Introduction 

The simulation of realistic models of cell membranes is of the utmost importance, 
if not for immediate medical use, certainly for the development of new conceptual 
and practical tools in the application of Molecular Dynamics (MD) methods to 
these and similarly complex systems. A lot of research activity has gone in this 
direction (see for instance |jl], |2|, |3|, ^ |5|, |6| and references therein) , but we are still 
far from having a complete understanding of the dynamic and thermodynamic 
properties of this kind of systems, not to mention the problem of simulating the 
formation and the dynamics of pores and ion channels (see, however, the recent 
very interesting work of ref. 0]). 

To predict the structure of molecular layers in different physico-chemical con- 
ditions by computer simulations, several modeling strategies have been proposed 
and investigated. Among them we would like to mention the following. 

1) Atomistic detailed approaches in which bilayer constituent molecules, water 
and possible counter-ions are modeled explicitly in full detail. The interaction be- 
tween electrostatic point charges are evaluated through Ewald related summation 
techniques, while the bias due to the use of periodic boundary conditions is re- 
duced by performing simulations in extended ensembles. Recent examples of this 
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kind of investigations can be found in refs. [Q, ^, 

2) Simplified models of structures formed by surfactants have been constructed 
that retain some of the essential interactions between the constituent particles, like 
the Lennard-Jones (LJ) forces among the hydrophobic tails and between solvent 
and hydrophilic molecular heads. On top of them an ad hoc repulsive soft core 
describing the interaction between hydrophobic and hydrophilic particles is added 
to allow self-assembling of micelles and bilayers 

3) Models in which only a single molecule or chain is modeled explicitly in 
terms of a realistic all-atom description, while the effect of the remainder of the 
system is parametrized by a mean-field deterministic force plus appropriate random 
forces. Models of this type are successfully treated by Monte Carlo simulations 
and/or stochastic techniques. In this context we wish to recall the fairly good 
agreement with experiments obtained in the evaluation of the NMR deuterium 
order parameters in the case of lipid bilayers, when Brownian Dynamics is used 
for their simulations |llO| . 

4) In the field of liquid crystals, where certain properties of layers often can be 
explored only through computer simulations, even more simplified force- fields, like 
that proposed by Gay and Berne (GB hereafter) with the successive modifi- 
cations of refs. and [13|, are used. Despite their simplicity these models allow 
the exploration and the study of various phases of the system, including smectic 
phases characterized by a number of different tilt angles. 

Comparing different models and trying to assess the effects of the various 
contributions that are included or neglected in different instances, is not a simple 
task. An example of this situation is the comparison between the GB potential |11] 
and its closest LJ representation. The GB potential was originally devised so as 
to mimic the interaction between two identical linear arrays composed by four 
LJ sites each. The properties of a fluid made of GB particles turn out to be 
rather different from its LJ counterpart: the GB fluid has a phase diagram with 
a well characterized nematic-isotropic transition temperature, while apparently 
the LJ counterpart does not show any phase transition of this kind when only 
the temperature is varied. This is a case where the use of a simplified potential, 
far from leading to an impoverishment of the model, turns out to give rise to an 
unexpected and very interesting thermodynamic behaviour. 

In the case of lipid molecules the comparison among different modeling strate- 
gies are much more difficult because of the occurrence of different types of in- 
teraction sites (head, tail and solvent) and because of the special importance of 
the internal molecular degrees of freedom, that may greatly influence the inter- 
molecular interactions, as well as the related hydrophobic tail packing and bilayer 
stability properties. 

The purpose of this paper is to perform MD simulations of a simplified, but still 
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sufficiently rich model of molecular bilayer in vacuum. The bilayer constituents 
are branched molecules composed by LJ sites with the topology and the inter- 
nal flexibility of the two fatty acid tails and glycerol group of the Dimyristoyl- 
phosphatidylcholine (DMPC hereafter) molecule. In the absence of water, the 
charged part of the DMPC head has been replaced by a methyl group (see Fig. |l|). 




Figure 1: The chemical structure of the DMPC molecule and the substitution 
which transforms it into what we have called DMMG. Hydrogens are not 
drawn. Carbons in the CH„ groups are labeled progressively in the standard 
way. 

With respect to an all-atom description of DMPC bilayers [Q, we retain in 
the model only the LJ interactions between the most hydrophobic portions of 
the lipid molecules together their intra-molecular flexibility. These simplifications 
allow us to single out and investigate the role of repulsive-dispersive interactions in 
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the structure of the bilayer, without the comphcations arising when electrostatic 
interactions among heads and solvent are present, which substantially limit the 
length of the simulated trajectories. 

Owing to these circumstances we are able to generate very long trajectories 
(up to ~10 ns) of sufficiently large systems, which allow us not only to expose 
and control certain effects related to the slow equilibration of the system (which 
would be important sources of systematic errors in the analysis of simulations not 
sufficiently longer than 1 ns), but also to accurately check the physical consistency 
of the parameters of the employed intra-molecular potential. 

We notice that in the simplified approach we present in this paper the bilayer 
components are not really amphiphilic. This is not, however, a completely silly 
approximation. We have checked, in fact, that the system configurations, occurring 
at the lowest density we have explored, in which it happens that some of the 
molecule tails are turned up-side-down, have approximately the same potential 
energy as the physical situation in which the tails of the molecules of the two 
layers are in contact. We find instead that at the surface density of both the 
crystal and gel phases the energy barrier which should be overcome by a molecule 
to rotate up-side-down by 180° degrees is sufficiently high to prevent this event 
to occur in the whole range of temperatures we have explored. This observation 
makes us confident that the predictions we have for the structural properties of 
the system in these two phases are fully reliable. As we said the situation, is 
completely different at our lowest density, which corresponds to that of the liquid- 
crystal phase. Evidently in this case the potential barrier is not sufficiently high, 
as we witness the crossing of it by some of the molecules. Despite this unphysical 
feature, we remark that our bilayer model remains stable at all the values of 
temperature and density we have studied (even at the lowest one), in the sense 
that, thanks to the strong attractive interactions between the two layers, none of 
the molecules leaves the layers. 

We also notice that at low surface density the potential energy of the whole 
system is presumably smaller than it would be if the repulsion between solvent and 
tails were introduced. In models in which solvatation effects are included |^ ^ low 
density configurations are essentially inaccessible, precisely because of a substantial 
solvent-tail repulsion. 

Given the crudeness of our bilayer model, we do not expect to be able to 
accurately predict thermodynamic parameters, like the dependence of the surface 
density or the magnitude of the stress tensor |^, |^ as functions of the temperature, 
mainly because of the lack of the interactions that govern the structure of the 
solvent interface and the absence of terms in the force-field that would prevent the 
tendency to isotropization at much too low surface densities we have mentioned 
above. Therefore, in order to study the physics of the different phases of our bilayer 
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model, we decided to fix the surface density at the experimentally estimated values 
(A'^yE'-simulations), rather than trying to work in an extended ensemble at fixed 
(lateral) pressure. 

As for the role of water, we tentatively assume that, once the bilayer is formed 
in the polar fluid, its structural and thermodynamic behaviour will be mainly 
governed by the packing and flexibility properties of the hydrocarbon tails |14]. 
Most of the effects of water on the detailed physico-chemical behaviour of the 
membrane are indeed of secondary importance for the kind of questions we are 
addressing here. In the picture we have in mind, in fact, water is only important in 
that it behaves as a sticky medium for the lipid heads, which modulates the surface 
density by more or less penetrating into the bilayer |15|. Thus water penetration 
is held responsible for the difference one sees in the increase of disorder in the head 
region, as compared to what happens to the tails, only to the extent it affects the 
surface density of the bilayer. The different behaviour between heads and tails 
may in turn be at the origin of the great sensitivity of the structure of the phase 
diagram to molecular composition and tail length. 

Even in the simplified model we are considering here, the known crystal struc- 
ture of the glycerol group and fatty acid tails of DMPC [16| is well accommodated 
in a local minimum of the force-field we are using, at the appropriate surface den- 
sity. The comparison of our results with those of more refined models can be of 
great help in understanding the relevance of more realistic modelizations of the 
molecular components of the system and/or the role of water. 

The MD simulations we have performed extend over a rather fine grid in the 
temperature-density plane. Thanks to this resolution and to the large statistics 
we have collected, we are able to draw interesting conclusions on the role of LJ 
interactions and internal flexibility on the stability and the structural properties 
of the crystal and gel phases of our bilayer model. 

The paper is organized as follows. In section 2 we describe the model of bilayer 
we have constructed and investigated. In section 3 we discuss the details of our 
simulation strategy and in section 4 we report the most significant results we have 
obtained with a special emphasis on the information we have collected concerning 
structure and properties of the crystal and gel phases of the system. Conclusions 
can be found in section 5. 



2 Simulation set-up 

Schematically, a cell membrane can be described as an almost spherical bilayer 
consisting of phospho-lipid molecules which separates the interior of the cell from 
the external world. Phospho-lipids are molecules composed by a hydrophilic head 
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and one or two hydrophobic tails. These pecuhar hydropathicity properties lead 
to the well known two-sheet structure of the membrane, in which the hydrophilic 
heads are in contact with water, present both outside and inside the cell, while the 
hydrophobic tails are more or less tail-to-tail pair-wise aligned. 

It appears experimentally that an important parameter governing the reaction 
rate of many biological processes, taking place inside or in the near vicinity of a 
membrane, is its "permeability" [14|. Prom this point of view a membrane can be 
thought as a system made out of two layers of a smectic fluid, with a permeability 
which depends "critically" on temperature, density, the detailed chemical compo- 
sition of the constituent phospho- lipids, the concentration of chemicals possibly 
dispersed in the membrane itself or in the solvent, etc. 

It is clear that a detailed simulation of the dynamics of the membrane of a 
living cell is just impossible and one has to resort to a number of simplifications. 
A fortunate circumstance in this respect is that it appears experimentally that 
the nature and the location of the phase transitions, which control a number of 
important physico-chemical properties of the membrane, are essentially related to 
the bulk ordering properties of the hydrophobic tails [14]. Thus a first step in the 
direction of simulating a realistic system, and the one which has been also largely 
followed in the literature |2|, ^, ^, is to take a two sheet system, each consisting 
of the largest possible number of lipid molecules (compatibly with the available 
computer power), in presence or even in absence of water, and proceed to study 
how the relevant order parameters behave as functions of temperature and surface 
density and what are the structural properties of the system in the different regions 
of its phase diagram. 

Lacking any realistic theoretical modeling of mesoscopic systems, as membrane 
or other molecular aggregates of biological relevance, one must resort to numerical 
methods, like MD or Monte Carlo and stochastic simulations of various kinds, 
in order to study the dynamic and thermodynamic behaviour of such complex 
systems. 

MD provides us with a microscopic, most often, classical |^ description of the 
system, which consists in following its time evolution by solving numerically the 
Newton equations of motion of its elementary constituents. 

The general mathematical setting of MD for the simulation of diffusive systems 
is well known [18] and we will not dwell on it here. For a general presentation of 
the method with a discussion of some of the tricks that have been developed to 
efficiently implement MD codes on the existing most powerful parallel platforms 
in the case of (NVE) simulations, we refer the reader to [p^]. 

We have prepared two chemically identical samples of largely different size 



^Car and Parrinello ]|r^ have proposed a consistent quantum mechanical generalization 
of the MD approach. 
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of our bilayer model, consisting of 2 x 32 and 2 x 256 molecules, hereafter de- 
noted for short BM{S) and BM{L) ("S"' and "L" stand for "small" and "large", 
respectively). Since, as we said in the Introduction, we do not have explicitly 
water in our simulations, we will neglect electrostatic interactions altogether. For 
consistency we have taken, as a model for the constituent molecules, an artifi- 
cial variant of DMPC, christened by us DMMG {Dimyristoyl-methyl- glycerol), in 
which the charged part of the DMPC head (defined as the part of the molecule 
above the first atom-a P atom-bound to the glycerol oxygen) has been replaced 
by a single CH3 group. The detailed chemical structure of the Y-shaped DMPC 
molecule is shown in Fig. where we have also indicated the modification we 
have made on its structure in order to construct the artificial variant of it we 
will be dealing with in this paper. The elementary constituents of the DMMG 
molecule will be schematically represented by LJ sites centered on the CH„ groups 
(n = 1,2,3), thereby using the so-called "united atom" approximation. In this 
way each DMMG molecule will consist of 37 elementary units (for short simply 
"atoms" in the following). Intra- molecular and fully flexible inter-molecular inter- 
actions are described by making reference to the OPLS force field pOl] with the 



modifications introduced in ref. |21]. According to the general OPLS philosophy, 
we have reduced the strength of the 1-4 LJ potential in each dihedral angle by a 
factor 8. The inter-molecular potential was linearly switched off between 1.0 and 
1.1 nm. As we already said, no electric charges are attributed to atoms. 

We have assumed the crystallographic density of our artificial system, po, to 
be the same as that of DMPC [16|. Experimentally the DMPC crystallographic 
surface density corresponds to an area of 0.396 nm^ per molecule in the plane of 
the bilayer (the x-y plane) . Despite the fact that we have drastically simplified the 
structure of the head of the constituent molecules, we decided to keep the same 
number of molecules per unit area as in the case of a DMPC bilayer, in order to be 
able to compare the properties of our model with those of a somehow related real 
system. At crystallographic density the surface of the system BM{S) spans an 
area of 3.56 x 3.56 nm^, that of the system BM{L) an area 8 times larger, equal 
to 14.24 X 7.12 nm2. 

The separation of the two layers in the z direction is fixed at the beginning 
of each simulation by setting the distance between the branching points (the C13 
atoms in the picture of Fig. |l]) of two opposite molecules, equal to 4.4 nm. 

In order to build the initial bilayer configuration we have started by generating 
a molecule with the dihedral angles of molecule A in Table 1 of ref [16|. Tail 7^1 
of this molecule ^ was taken to lie in the z-y plane, parallel to the z axis, with 
the head pointing in the positive z direction. Notice that tail #1 shows a kink 



^To be definite we call tail #1 the shortest of the two DMMG tails (the 014-C31 tail 
in Fig. P and tail #2 the longest one (the C15-C46 tail). 
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in correspondence to the dihedral angle O14-C17-C19-C20, making it parallel to 
tail #2 after C19. The basic molecule of the lower layer is generated by performing 
a rigid rotation of the first molecule by 180° around the y axis, holding fix the 
C13 atom. The two C13 atoms are then displaced by 4.4 nm, by translating one 
molecule with respect to the other one along the z axis. The fundamental cell of 
the initial crystal is an orthorhombic prism with the square basis in the x-y plane 
and the longer edge along the z axis. Two copies of this pair of molecules have 
been placed in the cell with the C13 atoms of the lowest molecules located one 
in a vertex and the other in the center of the x~y face, respectively. Finally the 
fundamental cell is replicated in the x and y directions the necessary number of 
times {i.e. 4x4 and 16x8 times, respectively). 

As usual, to limit finite-size effects, both BM{S) and BM{L) samples are 
taken to be periodic in the x and y directions. For simplicity of computation, 
periodicity is also imposed in the z direction, but adjacent copies of the bilayer 
are separated by about 100 nm, a distance so large that there are no possible 
interactions among different copies of the system. 

Two different MD codes, expressly developed by our group, were employed for 
the simulations of the two systems. Both codes make use of an MTS integration 
algorithm |22 , 23| with a long time step of 5 fs. A short time step 10 times smaller is 
introduced to integrate stretching and bending contributions. The package used for 
the BM{S) system (2 x 32 molecules) was run on a Digital 500 a-station. The code 
produces 1 ps of simulated dynamics every ~ 100 s of CPU-time (corresponding 
to 0.5 s of CPU time for every long time-step for a system of 2368 atoms). For 
the larger BM[L) system (2 x 256 molecules) we have employed a parallel version 
of the previous code. This code evolved from the one we used in |]l9| to test 
the level of the performances offered by parallel platforms in MD simulations of 
diffusive systems. Given the large request of computing power necessary to perform 
simulations of such a big system, use was made of the parallel platform called Torre, 
the largest (512 nodes, 25 Giga-flops) of the today available APE computers |24]. 
On this platform producing 1 ps of simulated dynamics costs ~ 2000 s of CPU- 
time (which corresponds to 2 s of CPU time for every long time-step for a system 
of 18944 atoms). 

In Tables || and § we report for the two samples {BM{S) and BM{L)) we have 
constructed the set of values of surface density and temperature (T in Kelvin), at 
which simulations have been carried out. The surface density is expressed in units 
of the crystallographic density po through the formula p = /i- po. The numbers in 
the entries of the Tables represent the length of the corresponding trajectory. 

At every density value, the sample BM{S) has been brought at the required 
temperature by a simple velocity-rescaling algorithm, which was run for more 
than 100 ps in the worst case. The surface density of the bilayer was fixed at 
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Table 1: The set of densities (p = ijl ■ Po, Po = crystallographic density) and 
equilibration temperatures (T) at which we have run the simulations of the 
system BM{S). The numbers in the entries of the Table represent the length 
of each trajectory in ns, not counting the initial velocity-rescaled steps. 



T{K)// p 
150 
225 
250 
275 
300 
325 
350 



1.05 1.0 0.95 0.91 0.87 0.83 0.79 0.76 0.69 



0.9 
0.9 
0.9 
5.1 
0.9 
0.9 
0.9 



5.4 
5.1 
5.1 
5.4 
5.1 
5.4 
5.1 



0.9 
0.9 
0.9 
0.9 
0.9 
0.9 
0.9 



0.9 
0.9 
0.9 
0.9 
0.9 
0.9 
0.9 



0.9 
0.9 
0.9 
0.9 
0.9 
0.9 
0.9 



5.4 
0.9 
0.9 
5.4 
0.9 
9.9 
5.4 



0.9 
0.9 
0.9 
0.9 
0.9 
9.9 
0.9 



0.9 
0.9 
0.9 
0.9 
0.9 
0.9 
0.9 



5.4 
0.9 
0.9 
5.4 
0.9 
5.4 
0.9 



the chosen value at the beginning of each simulation by uniformly rescaling the x 
and y coordinates of the center of mass of each molecule by the factor ^yji. The 
initial configuration of the larger system, BM{L), was obtained by replicating the 
last available configuration of BM{S) with the desired values of temperature and 
density, and then equilibrating the whole system for further 40 ps, again using 
a velocity-rescaling algorithm. Throughout this paper we will indicate for short 
by T the equilibration temperature, i.e. the temperature held fixed during the 
velocity-rescaling steps, and by Tg the real temperature at which the simulation 
was actually run. Of course the two temperatures will be somewhat different from 
each other. For most of our considerations in this paper this difference will be of 
little importance. 



Table 2: Same as in Table 1 for the system BM{L). Trajectory lengths are 
in ps. 



r(K)// p 


1.05 


1.0 


0.95 


0.87 


0.83 


0.79 


0.69 


150 


• 


140 


• 


• 


140 


• 


140 


225 


40 


60 


40 


40 


• 


40 


• 


275 


40 


140 


60 


40 


140 


40 


140 


300 


• 


60 


• 


• 


• 


• 


• 


325 


40 


140 


40 


40 


140 


40 


140 


350 


40 


60 


40 


40 


• 


40 


• 



10 



Few comments are in order here. 

1) The simulations of the smah system BM{S) are particularly long (they 
range from a minimum of 0.9 ns up to a maximum of 9.9 ns, not counting the 
initial velocity-rescaling equilibration steps) and span a grid of 7 x 9 = 63 points 
in the temperature-density plane (Table |l]). Notice that the trajectories employed 
in the analysis presented in the next sections are all longer than 5 ns. 

2) BM{L) simulations cover a similarly large region of temperatures and den- 
sities with a total of 28 points. Trajectories are, however, much shorter (Table ^). 



3 Simulation strategy 



One of the reasons to perform Molecular Dynamics simulations of such complex 
systems, like cell membranes, is that many biological processes are still waiting 
for an explanation at molecular level. Not many of the experimental techniques 
presently available may reach this goal at the resolution required and under really 
well controlled physico-chemical conditions. On the other hand, it is obvious that 
in order to be reasonably confident about the validity of the model used in a 
simulation, the structural properties of the model must be carefully compared to 
those of the "real" system so as to check the quality of the correspondence between 
simulation results and actual experimental data. 

In the case we are investigating, one of the most interesting macroscopic fea- 
tures characterizing the biological behaviour of the system is the structure of its 
phase diagram. The typical phase diagram of a phospho-lipid bilayer is repro- 
duced in Fig. |. The fi gure is taken from ref. ||2^ and refers to a Dimyristoyl- 



lecitin bilayer. As it is also seen from the figure, many phases have been identified 
experimentally: the most clearly separated ones are those called L, Lp and L^, cor- 
responding to the crystal, gel and liquid-crystal phase of the system, respectively. 
It is important to note here for the future comparison with simulation data that 
the horizontal axis in Fig. ^ is the percentage of water per phospho-lipid molecule 
and not directly the surface density. As we have argued before, however, to first 
approximation these two quantities are proportional to each other, as an increase 
of water concentration has mainly the effect of "diluting" the bi-dimensional fluid 
formed by the molecular heads. 

The nature of the structural changes associated to the different phase transi- 
tions has been mainly investigated in NMR experiments. NMR data |2^ indicate 
that, starting from the ordered crystal phase, the transition to the gel phase corre- 
sponds to an order-disorder transition of the heads, in which the tails retain their 
order, but undergo a collective "tilt" with respect to the plane of the bilayer [2 5( 1. 
The transition to the liquid-crystal phase is related to a further order-disorder 
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Figure 2: The typical phase diagram of a phospho-hpid bilayer. The figure 
is taken from ref. |^ and refers to a Dimyristoyl-lecitin bilayer. 
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transition, which this time involves also the tails. 

The experiments that put in evidence the existence of these different phases 
are generally performed on liposomes (sort of artificial small cells) in solution or 
on the so-called Langmuir-Blodgett multi-layers and are commonly performed at 
constant pressure. In these conditions the surface density of phospho-lipid heads 
changes with temperature. On the simulation side the requirement of holding the 
pressure constant can be fulfilled by performing MD simulations in the (NpT) en- 
semble. Simulations of this type are rather delicate as they require sophisticated 
strategies to force the system to move in the phase space with the appropriate sta- 
tistical weight and to tame otherwise annoying numerical instabilities. It should be 
immediately said, however, that today a lot of skill has been developed to properly 
deal with all these problems [^, |^, ^. Simulations in the (NVE) [micro- canonical) 
ensemble are definitely much simpler and numerical stability is easily achieved. 
In the micro-canonical ensemble pressure and temperature are not independent 
variables and can in principle be determined once the values of the density, N/V , 
and total energy, E, are given. 

Our strategy was to make a systematic study of the structural properties of 
the system by performing MD simulations at many different values of the {NVE) 
parameters in such a way as to explore quite a large set of values of temperatures 
and densities. The high "resolution" of our exploration has allowed us to draw 
interesting conclusions about the properties of the crystal and gel phase of the 
system. Results obtained for BM{S) are confirmed by the analysis of the data 
we collected for the much larger system, BM{L). As we shall see, BM{L) data 
are affected by almost identical statistical errors as the BM{S) data, though the 
available trajectories are much shorter. 

3.1 Extracting physical information from simulations 

Generally speaking, in order to show the existence of a phase transition one must 
be able to define an appropriate set of "order parameters" which could be used 
to describe the structure of the phase diagram of the system. The degree of or- 
der/disorder of the system can then be monitored through functions and parame- 
ters characterizing the translational and /or the orientational distribution functions 
of the constituent molecules. As a general reference for this kind of analysis, we 
shall follow in this paper the approach used in the discussion of liquid-crystal 
simulations in ref. p7| . 

As a first step in the study of the order properties of the system, we define the 
average direction of the constituent molecules. This quantity is represented by the 
1-st rank order parameter 

= iiu)f) (1) 
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where, for elongate molecules, like the ones we are dealing with, the unit vector u is 
naturally taken as the principal axis of the inertial tensor of the molecule [i.e. the 
normalized eigenvector of the inertial tensor with the lowest eigenvalue) . The dou- 
ble bracket in eq. ^ means average over both the Nmoi (identical) molecules of 
the system and over the N^onf configurations collected in the simulation. The 
superscript (^=1,2) means that the average is restricted to the ^-th layer (con- 
ventionally we decided to call layer #1 the upper layer and layer ^2 the lower 
one). The z component of P, which is usually called Pi, will be used to monitor 
the average orientation of the molecules. 

A further order parameter, which is also experimentally accessible (mainly in 
^H-NMR experiments [^), that can be used to parametrize the local orientational 
order of molecular segments, is 

5f) = ((i(3cos2A-l)))W (2) 

where the superscript has the same meaning as in eq. (|l]). In NMR experiments 
(3i is the angle between the z-th carbon-deuterium (Cj-D) bond in a tail and the 
external static magnetic field (which is usually taken to coincide with the normal 
to the bilayer surface). Si, often called "segmental order parameter", defines a 
good order parameter as it enjoys the following characteristic properties 

• When the Cj-D bonds are essentially aligned in a given direction with re- 
spect to the magnetic field, the system is in an ordered phase with Si ~ 
^ (3 cos^ — 1) / 0, where 13^'^^ is the most represented angular value in 
the distribution. In particular, if all the bonds have (3i = 0, then Si = 1. 

• When the distribution of cos /3j is flat [i.e. when all the values of cos /3j 
appear with equal probability), 5j = and the system is in a completely 
disordered phase. 

The crucial point that must be taken into account in trying to adapt the above 
definition to the situation one encounters in numerical simulations is that the order 
parameter we would like to define should allow us to distinguish in a clear way 
genuine order-disorder transitions from possibly present (and also very interesting) 
transitions between different kinds of ordered structures. The latter is precisely the 
situation that occurs experimentally in the crystal^gel transition, where, as we 
recalled above, the phase transition appears to be characterized by the following 
two structural changes: 1) heads undergo an order-disorder transition, 2) tails 
remain ordered, but the angle made by their average direction with the normal 
to the plane of the bilayer increases significantly (tilt). The problem with the 
definition (H) is that, in the situation in which the reference axis is taken to be the 
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same for all tail segments and lie parallel to the bilayer normal, an overall change 
in the direction of the tails with respect to it will actually affect the value of the 
order parameter (making it smaller), much in the same way a real increase of the 
disorder of the tails would do. 

To overcome this type of difficulties it was suggested in |^ that the way to 
appropriately monitor the local orientational order of atomic segments along the 
molecule is to use the formal expression given in eq. (^), but with /3j the Euler 
angles formed by the z axes of suitably chosen local reference frames (see below) 
with the average direction of the molecules in the layer to which the segments one 
is considering belong. Each local frame is defined by taking the z axis parallel to 
the segment connecting two non-contiguous carbon atoms along a tail. 

Let us now see how these general considerations are translated into precise 
formulae. Since we are interested in investigating 2-nd rank tensor properties of 
the system, as those that are extracted from NMR experiments, it is convenient to 



measure order with respect to the directions of "maximum 2-nd rank order" |18, 
^ . These directions are determined starting from the construction of the traceless 
tensor (whose definition closely parallels that of quadrupole moment of a charge 
distribution pi) 



QS = ((^(3^an^-(^a/3)))(') a,/3 = 1,2,3 (3) 



where -u is a unit vector related to the geometry of the molecules, b is the Kro- 
necker symbol and again one defines a Q-tensor for each one of the two layers. 
Eigenvalues and eigenvectors of Q reflect collective orientational properties of the 
system. For instance, if the phase is uniaxial (more mathematically, if the system 
is a collection of objects with D^a^h point space group symmetry), u should be 
taken as the direction of the principal axis of inertia of the molecules. Then Q 
will have one positive eigenvalue, say, A (usually called P2) and two identical neg- 
ative eigenvalues, —A/2. A biaxial phase is characterized by the fact that the two 
negative eigenvalues are not equal, while a more complicated orientational order 
is reflected in a non-trivial dependence of the eigenvectors of Q on the choice of u. 
In the uniaxial phase the (normalized) eigenvector corresponding to the highest 
eigenvalue. A, is called "director" and will be indicated by S-^'^ in the following. In 
our case c?^^ represents the average direction of alignment of the long molecular 
axes in layer L 

A precise notion of local order in each layer is introduced by comparing the 
orientation of a local rigid frame associated to any given molecular fragment with 
the orientation of a fixed reference frame, related to the collective order of the 
system, which we will identify with the (orthogonal) eigenvectors of Q^^^ Given 
three non-aligned atoms, indexed by zi, 12 and 23 the local frame associated to 
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each fragment, i (i = 111213), is constructed by taking the z axis as the vector 
joining the atom ii with the atom and the y-z plane as the plane where the 
three atoms lie. The azimuthal and polar angles (aj and Pi) of the director c?^) in 
these local frames can be used in turn to specify their orientation. 

It is important to remark that, in comparing values of the order parameter, 
Si, constructed through eq. (||) where the angles Pi are defined as described above, 
with numbers extracted from NMR experiments, the latter must be multiplied 
by a suitable numerical factor. This is related to the fact that in our definition 
of order parameter we monitor the orientation of the segment joining, say, the 
two non-contiguous Ci-Cj+2 atoms, while the two experimentally relevant Cj+i-D 
bonds, whose average direction is extracted from NMR data, lie in a plane almost 
perpendicular to it. Given the relative orientation of the different reference systems 
we have introduced, it can be shown, using Wigner rotation matrices, that in our 



geometrical conditions this factor is always very near to —2 [26|. 

To monitor the translational order / disorder properties it is customary to make 
reference to the pair distribution function. Given the strong anisotropy of the 
molecule arrangement, typical of our system, it is convenient to define the two pair 
distribution functions, and g±{r±), where ry and r± represent the length of 

the components of the vector joining the atoms of the pair in the direction parallel 
to the bilayer normal and in the plane perpendicular to it, respectively. As usual, 
and g± are normalized to the corresponding ideal gas distributions at the same 
density. We thus write 

^11 = = 7^ 

where h{r) is the number of pairs of atoms at (parallel, or respectively perpendic- 
ular) distance equal to r and h^{r) is the corresponding quantity for the ideal gas 
at the same density. It should be noted that we separately define a g± distribution 
for each one of the two layers by correspondingly limiting the calculation to the 
appropriate set of molecules. 



3.2 Data analysis 

For an ergodic system micro- canonical averages are computed as time averages of 
infinitely long trajectories. Such time averages are approximated in practice by 
averages over the largest possible number of independent configurations one can 
collect (given the available computing power) and representative of the infinite set 
of configurations of the system. As for the statistical error to be attributed to 
these ensemble averages, one should take the associated standard deviation. 

A look at the time evolution of the BM{S) system and, in particular, at the 
behaviour of the inter-molecular energy {Einter) shows that a plateau is smoothly 
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reached in most cases only after about 2.5-3 ns. In view of this rather slow equi- 
libration time we prudentially decided to exclude from the analysis presented in 
the next section the first 3 ns of each MD trajectory. 

The situation is much better for the BM{L) system for two reasons. One is 
that, as we explained before, the initial configuration of each BM{L) simulation 
was obtained by replicating the (already well equilibrated) last available configura- 
tion of the BM{S) system and further proceeding with an equilibration of the full 
system. The second is that the equilibration time of a larger system is naturally 
shorter. 

To decide how distant in time two successive configurations should be in order 
to consider them as uncorrelated, we studied the auto-correlation function of the 
potential energy. In the case of the BM{S) system, we found that, after removing 
from the analysis the configurations referring to the first 3 ns of each trajectory for 
the reasons explained above, the auto-correlation function dies away in a few ps at 
all values of density and temperature. A similar situation occurs in the case of the 
BM{L) system (though all the available configurations, produced after the initial 
velocity-rescaling equilibration, were kept in the analysis). Accordingly, we then 
decided to record a configuration every 1000 long time-step iterations {i.e. every 
5 ps) in all cases. 

For the small system this means that we will have at our disposal not less 
than 480 uncorrelated configurations (and sometimes much more, up to 1380, see 
Table |), that can be used in the forthcoming analysis. For the larger system 
BM{L), where much shorter trajectories have been collected, averages will be 
taken over a significantly smaller number of configurations. Despite this fact, 
thanks to the much larger size of the sample BM{L) compared to that of BM{S), 
statistical errors turn out to be of comparable magnitude. 

4 Results 

Out of the many surface densities we have explored we will show data corre- 
sponding to the three cases, fi = 1.0, 0.83 and 0.69, that roughly represent the 
experimental densities of the three most significant phases of the DMPC bilayer, 
i.e. crystal, gel and liquid-crystal, respectively. For each of these densities we will 
discuss results at three different temperatures (Tables || and|2|). The temperatures 
we have considered were chosen so as to have them below, in between and above 
the two most significant phase transitions (crystal^gel and gel^liquid-crystal) 
the DMPC bilayer undergoes (see Fig. ^). 

Despite the fact that we do not regard the results we obtained at the lowest 
density as physically interesting for the description of the expected liquid-crystal 
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Table 3: Inter-molecular {Einter) and intra-molecular (Eintra) energies at 
some selected temperature-density points for the BM{S) system. Numbers 
in the column labeled by Ts are the measured temperatures at which the 
corresponding simulation has actually run. Errors are in brackets. 



T(K) 




Ts (K) 


Einter (kJ/mol) 


Eintra (kJ/mol) 


150 


1.00 


157 (2) 


-207.5 


(0.6) 


84.8 (1.0) 


150 


0.83 


151 (2) 


-201.0 


(0.6) 


85.6 (0.9) 


150 


0.69 


160 (2) 


-195.9 


(0.6) 


96.5 (1.0) 


275 


1.00 


283 (4) 


-204.1 


(1.0) 


131.7 (1.6) 


275 


0.83 


281 (3) 


-196.7 


(1.0) 


132.3 (1.6) 


275 


0.69 


294 (4) 


-186.4 


(1.0) 


143.9 (1.6) 


325 


1.00 


338 (4) 


-199.7 


(1.2) 


154.6 (1.9) 


325 


0.83 


341 (4) 


-186.7 


(1.2) 


161.0 (1.9) 


325 


0.69 


343 (5) 


-185.0 


(1.7) 


175.8 (2.1) 



phase of a bilayer, we will nevertheless briefly discuss them in this section, mainly 
to contrast these data with what we get at the more reliable values of the density, 
/i=1.0 and /u=0.83. This discussion may also be very useful to identify where and 
how our model should be improved. 

4.1 The system BM{S) 

We shall start our analysis by discussing the data we get from the simulations of 
the system BM{S), for which we have collected very long trajectories (see Table |I|). 
We will use the results coming from the BM{L) system to refine and confirm the 
physical picture that emerges. 

4.1.1 Thermodynamic properties 

The density and temperature dependence of Einter and P2 should show marked 
variations in correspondence to first or second order phase transitions. In Table ^ 
a smooth and rather flat dependence of Einter on T and ^ is, instead, seen ^. 
This may be partly due to a still too coarse grid of points in the temperature and 
density plane and partly to the finite size of the simulated systems. 

■^In the third column of Table |^ we have reported the actual physical teniperature, Tg, 
of each simulation. 
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Table 4: Values of the 1^* rank order parameter, Pi, at some selected 
temperature-density points, separately for the upper (#1) and the lower 
(7^2) layer of the BM{S) system. Errors are in brackets. 



T(K) 


Pi(layer #1) Pi(layer #2) 


150 1.00 
150 0.83 
150 0.69 


-0.909 (0.004) 0.891 (0.004) 
-0.775 (0.003) 0.763 (0.004) 
-0.637 (0.005) 0.627 (0.004) 


275 1.00 
275 0.83 
275 0.69 


-0.923 (0.005) 0.913 (0.005) 
-0.788 (0.005) 0.782 (0.005) 
-0.63 (0.02) 0.58 (0.02) 


325 1.00 
325 0.83 
325 0.69 


-0.914 (0.007) 0.924 (0.005) 
-0.795 (0.007) 0.800 (0.007) 
-0.46 (0.03) 0.73 (0.02) 



At all values of density and temperature the molar heat capacity per atom, 
Cv, which was computed according to |3^, shows values in the range 3-3.5 (in 
units of the gas constant, R), which are typical of a molecular liquid. 

As for the parameter Pi (eq. (|^) ) , we see from Table |^ that it is rather sensitive 
to the average orientation of the molecules with respect to the z-axis (which, 
we remember, coincides with the bilayer normal) and, as expected, it steadily 
decreases with the density. 



Table 5: Same as in Table § for the 2"'^ rank order parameter, P2 



r(K) 


^^ 


P2(layer #1) 


P2(layer #2) 


150 


1.00 


0.986 


(0.002) 


0.996 


(0.001) 


150 


0.83 


0.979 


(0.004) 


0.962 


(0.003) 


150 


0.69 


0.961 


(0.004) 


0.949 


(0.005) 


275 


1.00 


0.996 


(0.001) 


0.995 


(0.001) 


275 


0.83 


0.995 


(0.001) 


0.995 


(0.001) 


275 


0.69 


0.95 


(0.02) 


0.80 


(0.02) 


325 


1.00 


0.996 


(0.001) 


0.995 


(0.001) 


325 


0.83 


0.978 


(0.008) 


0.990 


(0.002) 


325 


0.69 


0.28 


(0.05) 


0.84 


(0.02) 



AVe remark that the values of the parameters Pi , P2 and dz (see Tables |, 
^ and |6|) for the two layers are consistent with each other within three (most 
often, two) standard deviations for the two highest values of fi {i.e. /u=1.0 and 
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0.83). There is instead a statistically significaiit difference, especially at the highest 
temperatures, if one compares values referring to the upper and lower layer at 
the lowest density. This is due to rare events in which some of the molecules 
in the upper/lower layer rotate upward/downward relatively to the bilayer plane. 
Unfortunately the characteristic times of these motions are too long compared to 
the length of the trajectories we have produced to be in position of adequately 
sampling this dynamics. It is then not surprising to find that, within the time 
window of our simulations, the two layers do not appear statistically identical. 



Table 6: Same as in Table ^ for the director z component, dz. 



T(K) /. 


4 (layer #1) 4 (layer #2) 


150 1.00 
150 0.83 
150 0.69 


0.913 (0.004) 0.892 (0.004) 
0.781 (0.003) 0.773 (0.004) 
0.646 (0.005) 0.639 (0.004) 


275 1.00 
275 0.83 
275 0.69 


0.924 (0.005) 0.914 (0.005) 
0.790 (0.005) 0.783 (0.005) 
0.67 (0.02) 0.1 (0.7) 


325 1.00 
325 0.83 
325 0.69 


0.915 (0.007) 0.926 (0.006) 
0.801 (0.007) 0.803 (0.007) 
0.6 (0.2) 0.83 (0.01) 



Although we do not see any sign of instability of our bilayer model at anyone 
of the values of density and temperature we have explored (we have performed 
simulations at densities as low as /x=0.69 and temperatures as high as T=350 K) 
or after very long simulation times (some of the trajectories are as long as ~10 ns), 
we will not consider data referring to /i=0.69 as physically significant, in view of 
the undesired behaviour described above. 

All the data presented in this section are consistent with the description of 
the bilayer as an oriented molecular liquid, if one excludes the density ;U=0.69, 
where the system starts to show an increasing degree of isotropy. As for other 
structural parameters, particularly the pair distribution functions (see below), they 
show a high level of translational order that is progressively lost when the density 
decreases or the temperature increases. Finally, it is worth noticing that decreasing 
the density affects the orientational order in a stronger way than increasing the 
temperature, indirectly suggesting that ordering is more sensitive to hydration 
than to heating. 
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4.1.2 Structural properties 

We start by discussing some general structural properties of the bilayer at selected 
density-temperature points, that are suggested by a simple graphical inspection 
of the geometry of the system. In Figs. |^ and § we show two views ((a) and (b)) 
of two snapshots of the larger bilayer, BM{L), at T=325 K and at the two den- 
sities, fi=l and 0.83, respectively. The snapshots are obtained using the program 
MOLMOL |31|. In panel (a) of Figs. ^ and ^ we show the side view of the bilayer 



as seen from the direction in the x-y plane from which the translational order of 
the upper layer is most clearly visible. In panel (b) of Figs. ^ and |^ we show the 
top view of the same pictures, obtained after rotating the bilayer by 90° around 
an axis orthogonal to both the z axis and the direction of "best view" identified 
above. As a whole, the figures show that the main effect of decreasing the density 
is a progressive tilt in the average orientation of the molecules with respect to the 
bilayer normal. 

Further information on the translational order properties of the system can be 
gained by looking at Figs. |5|and|6|, where we have plotted at the two (density, tem- 
perature) values considered above, namely (a)=(l, 325 K) and (b)=(0.83, 325 K)), 
the projections on the x-y plane of the segment joining the two atoms of the pair 
C10-C13 (Fig. I) and C24-C37 (Fig. |), as they evolve in time. The two pairs 
we consider are taken to belong to the upper layer. They differ in their average 
distance from the middle plane lying in between the two layers: the C10-C13 pair 
can be used to describe the orientation of the head, as the two atoms belong to the 
top portion of the lipid molecule, while the C24-C37 pair can be used to monitor 
the orientation of the imaginary rectangle roughly containing the two molecular 
tails, because the C24 and C37 atoms belong to different tails and are located 
at more or less the same height from the C13 joint. A remarkable feature, quite 
evident in all the figures, is that molecules undergo fiuctuations in local positional 
energy minima that are spatially well ordered at high density and progressively 
less ordered as the density decreases. An other interesting question that can be 
studied by analyzing Figs. ^ and |6| is the effect of density on the degree of molecule 
packing: the range of positional fiuctuations at /i=0.83 is smaller than that at the 
higher density, fi=l, both for heads (Fig. |5|) and tails (Fig. ^). This means that 
the larger space the molecules have at their disposal, going from higher to lower 
densities, does not directly infiuence their actual translational freedom. The rea- 
son for this behaviour may be due to the fact that a decrease of the interaction 
potential among the molecules in each layer (that occurs because molecules are 
more separated) is at least partly compensated by the stronger interaction poten- 
tial between the two layers (which follows from the fact that an increase of the 
molecular tilt angle decreases the distance between the two layers). 
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Figure 4: Same as in Fig. ^ at /i=0.83 and r=325 K. 
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Figure 5: The projections on the x-y plane of the vectors joining the pair of 
atoms C10-C13, as seen in selected configurations of the BM{S) trajectory 
at T=325 K and /i=l (a), /i=0.83 (b). 
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(a) 



(b) 




Figure 6: Same as in Fig. ^ for the C24-C37 pair of atoms. 
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Finally, an inspection to the geometry of single molecules in different density- 
temperature conditions shows that at /i=l, T=150 K, the longer C15-C46 tail 
(tail #2) lies essentially parallel to the head, while the 014-C31 tail (tail #1) 
shows a kink in correspondence to the dihedral angle O14-C17-C19-C20, making 
it parallel to tail #2 after C19. This geometric arrangement, that is present in 
the initial crystallographic structure |16], becomes "inverted" at smaller densities 
and even at /x=l at the highest temperature, T=325 K, in the sense that the 
kink appears instead in the longer tail ^^2. In this situation the two tails happen 
to have almost the same effective length, if counted from the common branching 
point, C13. The related torsional changes occurring in the neighborhood of the 
branching region seems to be a very important ingredient in arranging the internal 
molecular structure so as to comply with the tilting and packing of molecules 
taking place in each layer. 



4.1.3 Pair distribution functions 

The information obtained through the graphical investigations we have just de- 
scribed can be quantified in a useful way with the help of the pair distribution 
functions and the order parameters defined in the previous section. 

In Fig. the distribution of CIO (panel (a)) and C24 pairs (panel (b)) at 
T=325 K and for the three densities, /U=l, 0.83 and 0.69, are compared. CIO is 
the group that replaces the charged part of the head of the DMPC molecule in our 
DMMG model, while C24 is a group located in the deep hydrophobic region. We 
see a significant difference in the translational structure between heads and tails, 
which is present at all densities: tails (panel (b)) look definitely more ordered than 
heads (panel (a)). At /x=l both tails and heads are ordered (see Fig. ^ and the 
solid line in Fig. |^ (b)): the system appears to be in a crystal phase. Of course 
ordering of the heads is progressively lost by increasing the temperature, as seen 
by comparing Fig. ^ (a) with Fig. ^. At /i=0.83 tails are still ordered but heads 
shows a structure which looks more liquid-like. This is very much indicative of a 
gel phase. Finally at the lowest density, /x=0.69, we see a strong disordering of the 
tails, owing to molecular upward/downward rotation, with a consequent tendency 
to isotropization. Looking at the shape of the dotted curve in Fig. (b) and at 
the Pi and P2 data at ;U=0.69 in Tables |^ and H, we are thus led to the conclusion 
that at this density the system finds itself in a configuration which is much too 
isotropic to be possibly interpreted as the liquid-crystal phase of a bilayer. It 
should be noticed, however, that tail disordering does not seem to affect the value 
of the inter-molecular energy, because, as we noticed above, the potential energy 
loss in each layer is compensated by a stronger interaction energy between the two 
layers. 
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Figure 7: The g± pair distribution function of CIO (a) and C24 (b) atoms for 
the system BM{S) at T=325 K and /x=l (sohd Hne), /i=0.83 (dashed hne), 
//=0.69 (dotted hne). 
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Figure 8: The g± pair distribution function of CIO atoms for the system 
BM{S) at T=150 K and (sohd hne), //=0.83 (dashed hne). 
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Figure 9: The g\\ pair distribution function of C13 atoms for the system 
BM{S) at T=325 K and ii=l (sohd hne), //=0.83 (dashed hne), /i=0.69 
(dotted hne). 
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There is a correlation among the decrease in density, the tilt of the molecules 
in the two layers and the decrease in the bilayer thickness. In Fig. ^ we plot 
the g'li pair distribution functions of C13 atoms, again for r=325 K at the three 
decreasing densities, /i=l (solid line), 0.83 (dashed line) and 0.69 (dotted line). 
The first thing that should be noted is the presence of very pronounced peaks at 
r|| ~ 0. These peaks come from the contributions of pairs of atoms belonging to 
the same layer. Concentrating our attention to the region ry > 0.5 nm, we see 
that the solid and dashed curves show narrow peaks, sticking out over a region 
of vanishing (^y, centered at ry ~ 3.8 nm and ry ~ 3.3 nm, respectively. They 
are indicative for the existence of two well separated layers. We remark that the 
first number is not too far from the value of the experimental bilayer thickness, 
which is ~ 4 nm at the same density and temperature. The shift in the location of 
this peak is nicely explained by the progressive molecule tail tilting, taking place 
by lowering the density, which tends to reduce the bilayer thickness. The overall 
shape of the dotted curve in Fig. ^ is, instead, completely different: first of all the 
rightmost peak is much broader than the corresponding peak seen in the other 
two curves and secondly there is no region where g\\ is zero: in fact, pretty soon, 
to the left of the rightmost peak (7|| starts to grow again. This behaviour reflects 
the fact that at this density the two layers partially overlap, as a consequence of 
the tendency to isotropization described above. Again, we notice that the effect 
of increasing the temperature on the bilayer thickness is negligible, if compared to 
what happens when we lower the density. For instance, increasing the temperature 
from 150 K to 325 K produces a thickness variation (which is maximal at ;U=0.83) 
of only 0.3 nm, while by decreasing the density from = 1 to = 0.83 and from 
/i = 0.83 to /i = 0.69, the peak in Fig. ^ gets shifted by 0.5 nm and 1 nm, 
respectively. 



4.1.4 Angular Distributions 

In Fig. |l^ we show the distribution function P{9) of the polar angle, 9, formed by 
the bilayer normal with the z axis of the local frame identified by the atoms C34- 
C40-C46 of tail ^^2. This direction, which, we recall, is chosen to be parallel to the 
vector joining the C34-C46 atoms, is practically aligned with the long molecular 
axis. We show the results for the lower layer only, as data for the other layer are 
essentially identical. The progressive tail tilting of the molecules with decreasing 
density is clearly visible. We also see that the peak of the distribution moves 
from 25° at /i=l (long-dashed line) to 40° at /i=0.83 (short-dashed line). At the 
lowest density (dotted line) the distribution becomes rather broad, but still peaked 
around 40°. The experimental tilt angle for similar systems in the gel phase, which 



has been measured to be about 30° [Q, is well consistent with these results. We 
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Figure 10: The distribution function at T=325 K of the angle, 9, formed by 
the segment C34-C46 of tail #2 with respect to the bilayer normal, for the 
BM{S) system at ^=1 (long-dashed hne), //=0.83 (short-dashed line) and 
//=0.69 (dotted line). 
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regard this agreement as a success of our model. 

A final structural investigation has been performed on the population of the 
different torsional states along the tails. Since the local torsional potential has 
three energy wells, corresponding to two gauche states ((/>=60° and 300°) with 
the same energy and one trans state ((/>=180°) with a different energy we can 
adequately characterize the situation by a parameter, usually called Pg, which 
represents the percentage of the total gauche population belonging to each torsion. 
We find at T=325 K the results reported in Table ^. In this computation we have 
averaged data over the two tails and the two layers. Values smaller than 1 % are 
not reported. 



Table 7: Percentage (%) of population of tail torsions at T=325 K for the 
BM{S) system at the various densities. Values smaller than 1% are not 
reported. 



Torsion 




H=0.83 


^=0.69 


(C12-C13-014-C17)+(C13-C15-016-C32) 


23 


25 


24 


(C13-014-C17-C19)+(C15-016-C32-C34) 


4 


4 


7 


(O14-C17-C19-C20)+(O16-C32-C34-C35) 


50 


50 


44 


(C17-C19-C20-C21)+(C32-C34-C35-C36) 


1 


10 


7 


(C19-C20-C21-C22)+(C34-C35-C36-C37) 




2 


6 


(C20-C21-C22-C23) + (C35-C36-C37-C38) 




2 


5 


(C21-C22-C23-C24) + (C36-C37-C38-C39) 




1 


5 


(C22-C23-C24-C25) + (C37-C38-C39-C40) 




1 


5 


(C23-C24-C25-C26) + (C38-C39-C40-C41) 




1 


4 


(C24-C25-C26-C27) + (C39-C40-C41-C42) 




1 


5 


(C25-C26-C27-C28) + (C40-C41-C42-C43) 




1 


5 


(C26-C27-C28-C29) + (C41-C42-C43-C44) 




1 


6 


(C27-C28-C29-C30) + (C42-C43-C44-C45) 


1 


4 


8 


(C28-C29-C30-C31) + (C43-C44-C45-C46) 


3 


6 


9 



At the two highest densities only the last two torsions (C27-C28-C29-C30, C28- 
C29-C30-C31) in tail #1 and the corresponding ones in tail #2 show a population 
larger than 1%. The high values, Pg=2b% and Pg=50%, that we find for the 
average percentage population in the torsions (C12-C13-014-C17) + (C13-C15- 

^Actually the energy of the trans state is lower than the energy in the gauche states in 
all torsions involving only aliphatic carbons. 
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016-C32) and (O14-C17-C19-C20) + (016-C32-C34-C35), respectively, are due 
to the distortion of the geometry in the "kinked" region close to the ester bond. 
These values are essentially the same for all densities and temperatures and are 
due to the phenomenon of "kink inversion" between the two tails, described before. 
It is worth noticing that the ^=0.83 results are consistent with those obtained for 
the DPPC in the gel phase after correcting for the fact that the DPPC 

hydrophobic tails are 2 methylene groups longer than the DMPC tails. 

The smooth increase of the gauche population along the tails, seen in the 
simulations performed at the lowest density, means that the rotational freedom of 
molecule tails does not help in making the tail flexibility in the terminal region 
larger than that in the core. This is at variance with the results of the first paper 
of ref. IQ, where a larger mobility of tail end segments, as compared to what is 
seen in the core region, was found in simulations of the liquid-crystal phase of 
DPPC. Values around 20 % for Pg for the last six torsions of the tails are, in fact, 
reported. One might argue that the little tail flexibility we find in our model might 
be attributed to the much too strong 1-4 LJ potential we have introduced in each 
torsion. These 1-4 forces mainly have the effect of adding a steric repulsion in the 
gauche conformations. We recall that within the framework of the OPLS force- 
field |2^, which we employed in the present paper, we have reduced, according 
to the general OPLS philosophy, the strength of the 1-4 LJ potential by a factor 
8. This may not have been enough in view of the too high tail rigidity we have 
found. However, neglecting 1-4 forces altogether gives results only partially in the 
desired direction. In fact, on the one hand it actually produces an increase of 
the gauche population by a large factor (in trial simulations performed on liquid 
butane this factor is found to be about 3), but on the other hand it leads to 
a more pronounced upward/downward molecular rotational freedom, which may 
make even more problematic for the system to unveil its liquid-crystal phase. 



4.1.5 Segmental Order Parameter 



In Fig. |ll| we show the behaviour of the order parameter, Si (eq. @ of sect. 3.1), as 
function of the position of the (Cj-Cj+2) segment along the tail ^ Panels (a) and 
(b) show the results for tail #1 and #2, respectively, averaged over the two layers. 
We focus on the behaviour of Si at r=325 K and different densities. We see that, 
correctly, the order parameter decreases with the density (the expectation for a 
completely isotropic fluid would be 5j=0). Notice that, not surprisingly, data at 
the lowest density are affected by errors that are definitely larger than those on 



^We observe that the notation (Ci-Ci+2) does not apply to the initial segments of the 
two tails, that are called C17-C20 and C32-C35 respectively, because in the enumeration 
of carbons the labels 18 and 33 are attributed to O atoms. 
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the data at /x=l and /i=0.83. 




Figure 11: The S order parameter for the system BM{S) at T=325 K as 
a function of the position along tail #1 (a) and tail #2 (b): stars are for 
/i=l, squares for /i=0.83 and triangles for /i=0.69. The first segment starts 
from the carbonyl carbon atom C17 in tail #1 and C32 in tail #2. The lines 
joining the points are only to guide the eye. 



4.2 The system BM{L) 

As summarized in Table we have performed extensive MD simulations also on the 
much larger system, BM{L), consisting of 2 layers of 256 DMMG molecules each. 
To compare with previous results we have analyzed the trajectories corresponding 
to the same combined values of surface density and temperature considered above. 
An immediate observation that emerges from the analysis of BM{L) MD data is 
that statistical errors are about the same as those one finds in the case of the system 
BM{S), despite the fact that substantially shorter trajectories were produced 
(compare Tables ^ and |2|). As we already noticed, this is a direct consequence of 
the increase in the size of the system and of the fact that the larger BM{L) system 
was built by assembling already equilibrated BM{S) blocks with the result that 
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no special needs for further long or sophisticated equilibration procedures have 
emerged. 

The general conclusions of the analysis of BM{L) data is that the scenario we 
have come up with in the previous sections is well confirmed by the simulations 
of this larger system, and we have good consistency between the two sets of data 
within errors. 
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Figure 12: The S order parameters at r=325 K and ^=1 as function of the 
position along tail #1 (a) and tail #2 (b) for the BM{L) system (squares) 
and for the BM{S) system (triangles). The lines joining the points are only 
to guide the eye. 

As an example of the quality of the results that we have obtained, we compare 



in Figs. |lj and |13| the behaviour of the order parameter S, for the two systems, 
BM{S) and BM{L) at T=325 K and /i=l and /x=0.83, separately for the two tails, 
averaged over the two layers. The BM{L) points are systematically lower than 
the BM{S) data, as expected from the lowering of the collective orientational 
order with the increasing size of the system. This effect is more relevant for 
segments close to the head and to the end of the tails. The slightly larger errors 
that affect the data associated to the first and the last tail segment are due to 
jumps among different torsional states that are not sufficiently well sampled in 
the time evolution at our disposal. In order to reduce errors associated to this 
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slow torsional dynamics, it would be important to increase the length of BM{L) 
simulations. Data at ;U=0.69 are affected by substantially larger errors and are not 
shown. 
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Figure 13: Same as in Fig. O at yU=0.83. 



The above indications about the behaviour of the degree of collective orienta- 
tional order with the size of the system are confirmed by a careful study of the 
order parameters, Pi and P2- For instance, at ^=0.83 and r=325 K one gets 
for the BM{L) system P2=0.964 it 0.003, a value somewhat lower than the one 
observed for the smaller BM{S) system, for which an average value P2=0.984 it 
0.002 is found. A lowering of the orientational order with the size of the system is 
actually expected from the observation that in a bigger system fluctuations with 
longer wave lengths are possible. 

Further agreement among BM{L) and BM{S) simulations emerges from the 
study of translational order properties. The g±{r) distribution functions computed 
from the BM{L) simulations show essentially the same structure and locations 
of the peaks that is seen in the BM{S) distributions. The only difference is 
a somewhat less pronounced peak structure in distributions referring to atoms 
belonging to the hydrophobic core. Given the strong similarity of BM[L) and 
BM{S) results, we do not report the former here. 

We would like to conclude this section with some general considerations about 
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the quest for MD simulations of larger and larger systems. Given the satisfactory 
level of consistency we find by comparing the MD data of the two samples, one 
might ask whether simulations of systems as large as BM{L) are really necessary 
to get the physics of the system right. Although the answer to this question 
strongly depends on the kind and the degree of detail of the information one 
wishes to obtain from the model, it is important to recall that the study of finite 
volume effects is a necessary step in any numerical simulation, because features 
and properties seen at small volumes may completely change when the volume is 
increased. Consequently it is absolutely mandatory to get an idea of the magnitude 
of finite volume corrections. Prom our experience it was only thanks to the data 
we had for the larger system that we could reliably draw useful conclusions on the 
structural properties of our model-membrane. 

Apart from these rather general considerations, there are in our opinion many 
specific problems where the fact of having a sufficiently large system is crucial in 
order to be in the position of extracting reliable results from numerical simulations. 
Just to make two examples we may mention: 

1) the problem of trying to accurately locate the critical lines of the phase 
diagram of the system. The obvious reason for the necessity of comparing systems 
of different sizes is that a structural change of the system can be interpreted as 
a true phase transition only if one is able to show that the transition becomes 
sharper and sharper as the volume increases. 

2) Another very interesting problem is the study of the dynamic and the for- 
mation of ion channels. In this case dealing with a sufficiently large system is 
required by the need of having the cross section of the channel sufficiently smaller 
than the surface of the bilayer itself in order to meet as accurately as possible 
the actual physical situation. We wish to mention in this context the pioneering 
investigation recently carried out in the beautiful work of ref. Q, where it has 
been shown that indeed channels have the possibility of forming and dynamically 
maintaining themselves stable. 



5 Conclusions 

We have performed a systematic MD study of the structural properties of a simple 
atomistic model of bilayer by analyzing the behaviour of a large number of physical 
quantities characterizing the system, as functions of density and temperature, and 
by comparing results obtained from simulations carried out on two systems of 
largely different size. 

We have found that our modelization reproduces in a remarkable way many 
known features of the gel phase of DMPC bilayers [25|, as well as the structural 
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changes that accompany the transition from the crystal to the gel phase. These 
findings have emerged quite clearly by analyzing the very long trajectories we have 
collected for the system, BM(S), and are nicely confirmed by the statistically 
equally accurate data we have for the much larger system, BM{L). 

Clearly a lot of work is still needed to further validate and improve our ap- 
proach, especially in the direction of investigating the role of water, whose effects 
were only indirectly accounted for through the modulation of the surface density 
of the system. The inability of the model to cope with the expected liquid-crystal 
phase of the system at low density is its major limitation. To prevent the extru- 
sion of the molecule tails from the planes of the bilayer at low density and high 
temperature, either some sort of geometric or dynamic constraint must be intro- 
duced, perhaps by modifying the intra-molecular tail flexibility, or more directly 
one should explicitly introduce water. 

Despite the limitations we have mentioned, from the quality of the results 
presented in this investigation, we conclude that the membrane model we have 
developed is capable of capturing most of the essential structural properties of the 
real system in the crystal and gel phase. Thanks to its simplicity, we are confi- 
dent that it will be possible to appropriately improve the model in the directions 
mentioned above, without affecting the key features that are at the basis of the 
stability of our bilayer, i.e. tail packing within the layers and adequately strong 
interactions between the two layers. 
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